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Abstract — A major effort in combustion research at the present time is devoted to the theoretical modeling 
of practical combustion systems. These include turbojet and ramjet air-breathing engines as well as 
ground-based gas-turbine power generating systems. The ability to use computational modeling exten- 
sively in designing these products not only saves time and money, but also helps designers meet the quite 
rigorous environmental standards that have been imposed on all combustion devices. The goal is to 
combine the very complex solution of the Navier-Stokes flow equations with realistic turbulence and 
heat-release models into a single computer code. Such a computational fluid-dynamic (CFD) code 
simulates the coupling of fluid mechanics with the chemistry of combustion to describe the practical 
devices. This paper will focus on the task of developing a simplified chemical model which can predict 
realistic heat-release rates as well as species composition profiles, and is also computationally rapid. We 
first discuss the mathematical techniques used to describe a complex, multistep fuel oxidation chemical 
reaction and develop a detailed mechanism for the process. We then show how this mechanism may be 
reduced and simplified to give an approximate model which adequately predicts heat release rates and a 
limited number of species composition profiles, but is computationally much faster than the original one. 
Only such a model can be incorporated into a CFD code without adding significantly to long computation 
times. Finally, we present some of the recent advances in the development of these simplified chemical 
mechanisms. 


INTRODUCTION 

One of the significant research thrusts in the field of 
combustion research today is the theoretical model- 
ing of practical combustion systems. These include 
turbojet and ramjet air-breathing engines as well as 
ground-based gas-turbine power generating systems. 
The ability to use computational modeling exten- 
sively in designing these devices not only saves time 
and money, but also helps us meet the quite rigorous 
environmental and safety standards which have been 
imposed on all combustion devices. With the develop- 
ment of high-speed parallel-processor computers sig- 
nificant advances have been made in the field of 
mathematical modeling of multidimensional turbu- 
lent, reactive flow processes. This modeling, of 
course, involves the numerical solution of the 
Navier-Stokes flow equations along with models for 
turbulence generation and heat-release in the flow. M 
These computational fluid-dynamic (CFD) codes 
generally require large amounts of computer time for 
a complete solution, even with today’s advanced 
machines. 

In practical combustion systems the process of 
heat-release is coupled with the flow phenomena, 
making a complicated modeling problem much more 
difficult. For many years combustion effects on the 
flow process were considered by using very simplified 
and unrealistic chemical reaction approximations in 
CFD codes. The combustion process itself is very 
complicated and requires a large size numerical- 
analysis computer code to be modeled accurately. 5,6 


Even with today’s computer technology, the linking 
of a fluid-dynamic code with a full combustion 
kinetics analysis code is not feasible because of the 
excessively long computational times which would 
result. Therefore, research in recent years has been 
directed towards developing simplified, computation- 
ally rapid, but still realistic models for chemical 
reaction and heat-release. 7 ' 9 Only this type of model 
can be incorporated into existing flow analysis codes 
without adding significantly to long computational 
times. 

In this paper we first discuss the equations which 
describe a complex, multireaction fuel oxidation pro- 
cess and the mathematical techniques used to solve 
them. We then describe a computer code, LSENS, 
which implements the most efficient solution method, 
and show its use in developing a detailed molecular 
reaction mechanism for the oxidation of a particular 
hydrocarbon fuel. 

The remainder of this paper focuses on methods of 
transforming a detailed mechanism into a form which 
adequately describes heat-release rates and a limited 
number of composition profiles, but is computation- 
ally faster than the original mechanism. There are 
several approaches to this simplification task. Among 
them are (1) use of pseudo steady-state approxi- 
mations, (2) use of partial equilibrium and rate-con- 
trolled constrained equilibrium techniques, and (3) 
development of a quasi-global chemical mechanism 
in which most or all of the molecular reactions are 
replaced by a smaller number of overall or global 
steps, whose rates are determined empirically by 
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fitting computed results to a fairly limited range of 
experimental laboratory measurements. A relatively 
new approach to simplifying chemical kinetic compu- 
tations is “solution mapping”, which attempts to 
replace the numerical solution of a set of differential 
equations by a group of algebraic equations. We will 
discuss these various approaches to the problem of 
computational simplification and present the current 
status of recent work in this field. 

The present paper is not meant to be a comprehen- 
sive discussion of all chemical mechanism reduction 
techniques, but only those which are applicable to the 
objective of predicting system responses in practical 
problems such as CFD modeling. We do not, there- 
fore, consider the matter of simplification by math- 
ematical or chemical “lumping” of reactions. These 
methods are described by Frenklach 8 and are appli- 
cable to exploratory problems such as determining 
the mechanism of soot formation in combustion 
systems. 


is a function of temperature only and is computed 
from fundamental thermodynamic relationships. 

Reaction rates and species formation rates 

For our constant volume closed (fixed mass) react- 
ing system the mole number of mole per unit mass of 
species i present will be called a,. The molar rate of 
change of species i per unit mass is written as by 
Bittker 6 


= t 0 ) = O/o (4) 

where NRS is the number of reacting species (which 
may be less than NS if some species are inert), p is 
the gas mixture density and IF, is the molar formation 
rate of species i per unit volume given by 


EQUATIONS FOR CHEMICAL REACTION OF A 
HOMOGENEOUS GAS MIXTURE 

Chemical reactions and rate equations 

We consider the reaction of a homogeneous gas 
mixture of NS species which participate in NR simul- 
taneous molecular reactions which can be written 
symbolically in the general form 10 

ns k NS 

7 = 1,2,..., NR. (1) 

i-i i=i 

Here v' ; is the number of moles (i.e. stoichiometric 
coefficient) of reactant species i in the / h reaction, v" 
is the number of moles of product species i in the / lh 
reaction and 5, represents the chemical formula of 
species i. If a species does not participate in a 
reaction, both v'j and v" are equal to zero. The arrows 
indicate that the reaction may be reversible and kj and 
k_j are the corresponding forward and reverse direc- 
tion rate coefficients for the reaction. Each kj is a 
function of temperature only and is usually given by 
the so-called modified Arrhenius expression 11 


NR 

w,= I CO,. (5) 

j - 1 

Here, co y is the net rate of formation of species i per 
unit volume due to reaction j, and is computed from 
the forward and reverse rates and stoichiometry of 
the reaction. For molecular reactions these rates obey 
the law of mass action. 12 This law says that a reaction 
rate is proportional to the product of the concen- 
trations of all reactants, each raised to a power equal 
to its stoichiometric coefficient in the reaction. 

For constant-density and either constant or as- 
signed temperature problems this set of ODEs com- 
pletely specifies the time evolution of the chemical 
system. The pressure, p, is computed from the ideal 
gas law 6 


Here, M„. is the mean molecular weight of the mixture 
which is related to the {cr,} by 


kj ~ Aj T n i exp^ — (2) 

Here A Jt n f and £) are constants, T is the Kelvin 
temperature and R is the universal gas constant. Each 
reaction may be either irreversible (i.e. proceed in the 
forward direction only) or reversible. In the latter 
case, the reverse rate coefficient, k_ jt is related to k } 
by the principle of microscopic reversibility 12 


where K cj is the concentration equilibrium constant 
for the y th reaction. For ideal gases this latter quantity 


= (7) 

i= 1 

If temperature is not assigned or density changes in 
a chemical reaction, differential equations are needed 
for these variables. Also, if reaction is in a flowing 
system, a differential equation for the rate of change 
of velocity would have to be obtained from the 
momentum conservation law. In this discussion we 
considered only a constant-volume static reaction in 
order to simplify the presentation of the theory. 
Therefore, only an additional temperature equation 
may be needed. 
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Temperature differential equation 

The variation of temperature, if it is not assigned, 
is obtained from the energy conservation law for a 
closed system 1 ’ 


(2)] must be known. Equally as important as the rate 
coefficient values is a data base of thermodynamic 
data which are used to compute the reverse rate 
coefficients [Eq. (3)]. 


d u + p dr = — d Q, (8) 


where u is the mixture internal energy per unit mass, 
Q is the heat-transfer rate per unit mass of mixture 
from the reacting gas, and v is the specific volume of 
the gas. We use the relation 6 

h = u + pv (9) 


and replace v by the reciprocal of density. Then, 
differentiation with respect to time gives 6 

d h 1 dp d Q 

— = — ( 10 ) 

dt p d t dt 

In these equations h is the enthalpy per unit mass of 
mixture. For an ideal gas h is given by 

NS 

h = I a,H n (11) 

1= l 

where //, is the molar enthalpy of species i and is a 
function of temperature only. Differentiation of Eq. 
(11) gives 


d h 
dt 


d T N « s 


dtr, 

d7 ’ 


( 12 ) 


Differential equations for sensitivity analysis 

When performing chemical kinetic computations 
an important question which must be answered is 
how changing the input parameters affects the com- 
puted results. This is especially true for the rate 
coefficients because their values are often not pre- 
cisely known. Research in the field of sensitivity 
analysis has developed methods of answering this 
question. 11 19 These methods are of two types: (1) 
local methods and (2) global methods. Local sensi- 
tivity analysis methods give the effect of relatively 
small changes of a parameter above and below its 
original value. We discuss in this paper only these 
methods, which have been highly developed. Global 
methods provide information about the effect of 
large-magnitude changes in any input parameter, but 
are more difficult to implement for practical use. We 
now give a brief description of the theory of local 
sensitivity analysis. 

The purpose of sensitivity analysis is to compute 
coefficients which tell the percentage change in any 
dependent variable caused by a given per cent change 
in an input parameter. We limit this discussion to 
changes in the three constants A n n, and E, in the rate 
coefficient expression [Eq. (2)] for any reaction. First, 
the previously derived ODEs for a reacting system are 
rewritten in a more general form. Equations (4) and 
(14) can be expressed as the set 


where 


NS 

c p = I cr,(C p ), (13) 

;= i 

is the specific heat of the mixture and (C r ), = dHJdt 
is the molar heat capacity of species ;, also a function 
of temperature only. Equation (12) can now be used, 
along with the perfect gas law, Eq. (6), to eliminate 
dh/dt from Eq. (10). The result of these substitutions, 
after simplification and rearrangement, is the follow- 


ing ODE for temperature 10 

dT 

~T = T-{ 
dt 

"ydo, H, dtr, Id Q 

,r, d 7 ~~ rt d7 - Rr~d7 

c r 1 


R M w 


T(t = t 0 ) = T 0 . (14) 

This equation is solved with the NRS equations of 
Eq. (4) for the NRS + 1 unknowns <r, , tr 2 , ... , o NRS , 
and T. In order to evaluate the dtr,/dt terms a 
chemical reaction mechanism must be formulated 
and expressions for the forward rate coefficients [Eq. 


^ =f{{yff{A l },{nff{E l }). (15) 


In this set of equations index / ranges from 1 to N, 
the number of dependent variables, and j goes from 
1 to NR , the number of chemical reactions. The first 
order sensitivity coefficient S, t is defined by 10 


.S',, = 


dqj 


k*.i 


(16) 


where rj , represents any of the Aj , n t and £} values. 
This coefficient is the rate of change of any computed 
y , when any rate parameter is changed. Since we are 
considering only a small change Ay in the parameter, 
the corresponding change in y\ can be approximated 
by multiplying the latter quantity by the sensitivity 
coefficient, S,j. 

A set of ODEs for these sensitivity coefficients can 
be obtained by first differentiating Eq. (15) with 
respect to rj, to get 


_£_/drA y, ffdy, 
8ti,\dt ) /“] cy t dtp 


Vk 

I A- j 


or], 


( 17 ) 
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The order of differentiation with respect to t and tjj 
may be interchanged and the definition of Eq. (16) 
used to obtain 


d Sjj 
d t 


y M 

/- 1 Sy, 



all 




(18) 


The initial conditions are S y (i = t n ) = 0.0, 
i = l,2,..., iV. Although the equation set [Eq. (15)] 
is quite nonlinear, the above equations for the {S u } 
are linear ODEs. They are, therefore, easier to solve. 
This solution can, in fact, be accomplished at the 
same time that the chemical modeling equations [Eq. 
(15)] are solved. 

Sensitivity coefficients are generally presented in 
normalized form so that meaningful comparisons can 
be made for widely different values of For any 
nonzero t]j this normalizing is easily done by using the 
relation 


iljdy, 8 In y t 
ij y, 5r\j 6 In r\j ’ 


(19) 


where <S, 7 > is the normalized sensitivity coefficient, 
or the per cent change in y, caused be a 1 % change 
in rfj only. This normalization is satisfactory for 
t]j = Aj, since this rate coefficient parameter is always 
nonzero. However, Eq. (19) cannot be used for the n t 
and Ej rate parameters, because these may be zero. 
Other normalization procedures have been devel- 
oped 10 and in all cases <S, 7 > represents the approxi- 
mate change in y, caused by a 1 % change in the actual 
rate coefficient kj, which is not always a 1% change 
in the rate parameter of interest. 


SOLUTION OF KINETICS AND SENSITIVITY ANALYSIS 
EQUATIONS 

Numerical solution of chemical kinetics ODEs 

The first attempts to solve numerically the ODEs 
modeling complex combustion reactions showed that 
classical methods such as the explicit Runge-Kutta 
and Adams methods were highly inefficient. The very 
small stepsizes these methods require are due to the 
well known property of “stiffness” shown by these 
and many other systems of nonlinear differential 
equations. 5,20,21 Many implicit numerical integration 
methods for solving such stiff equation systems have 
been presented over the last 40 years. The first to be 
applied to the solution of chemical rate equations was 
the backward differentiation method of Curtis and 
Hirschfelder. 22 Many other techniques have been 
described. 6,23 " 31 The efficiency and accuracy of several 
of these methods, as applied to combustion problems, 
have been examined by Radhakrishnan. 32 34 This 
work showed that the code LSODE 30 ' 31 is the most 
accurate and efficient algorithm now available for 
solving these stiff differential equation systems. 
LSODE is a package of numerical solvers for stiff and 
nonstiff ODE systems. For stiff problems the implicit 
backwards differentiation formula (BDF) method is 


used because of its property of stability for these 
problems. Thus, the step size is limited only by 
accuracy requirements and not by stability consider- 
ations. LSODE continuously varies the order of the 
BDF formula (up to a maximum of five) and the step 
size to obtain good computational efficiency. The 
code also allows a choice of several corrector methods 
in the predictor-corrector iteration procedure. 

We mention here two computer codes in common 
use for solving general chemical kinetics problems for 
any chemical system. The CHEMKIN code 35 allows 
the user to set up a chemical mechanism and then 
select from a library of subroutines which define the 
differential equation system to be solved. This code 
can be used with any numerical integrator, but 
LSODE is recommended and used by the authors. 
The NASA Lewis Research Center general kinetics 
code LSENS 10,36 also permits the use of any chemical 
system. However, it has the differential and algebraic 
equations for many different reaction models and the 
LSODE integrator built into it. The user does not 
select FORTRAN routines at execution time. The 
input data file allows many different problems and 
options to be selected quite simply. This code will be 
described in more detail in the following sections. 

Numerical solution of sensitivity analysis ODEs 

The most direct method of determining the effect 
of variations in rate coefficients on solution results is 
to change any rate parameter from its original value 
and repeat the calculation with the modeling code. 
This “brute force” approach is not only time consum- 
ing but also becomes very computationally expensive 
when the number of rate parameters is large. For this 
reason several methods have been developed for 
solving the system of equations, Eq. (18) for the { S,j } 
rather than computing them from the brute-force 
variations. The solution is obtained at the same time 
that the chemical kinetics equations are solved. Some 
of the methods that have been used are the direct 
method (DM), 37 the Green’s function method 
(GFM) 38 and its modifications 39,40 and the decoupled 
direct method (DDM). 41,42 

The DM solves the system of 2N simultaneous 
differential equations given by the system of 
equations, Eqs (15) and (18) in one operation for all 
the dependent variables and their sensitivity co- 
efficients. This technique was found to be rather 
inefficient and unstable and could not solve some stiff 
problems. 5 A variation of the direct method is the 
uncoupling of the model equations from the sensi- 
tivity coefficient ODEs. Equation (15) is solved first 
and the complete solution is stored for use in solving 
Eq. (18) separately. Although this technique should, 
in theory, be more efficient, it has shown no efficiency 
gains over the coupled direct method, and is also 
unstable for some stiff problems. 5 

The several versions of the GFM compute the 
Green’s function for the sensitivity equations and 
then obtain the {S, 7 } by integration over the Green’s 
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function. They also do the model solution first and 
then use the results to integrate the Green's function. 
Both the accuracy and efficiency of the GFM tech- 
niques are sensitive to the selection of error control 
parameters and the choice of printout stations for the 
answers. 41 

The DDM also solves the model equations separ- 
ately from the sensitivity coefficient equations. But 
the two solutions are done sequentially, i.e. step by 
step. The results of the model solution on any step are 
used to solve the sensitivity equations for that same 
step, so there is no need to store a large amount of 
information about the model solution. For solving 
the stiff problems of combustion chemistry the DDM 
has been shown by Dunker 41 to be more efficient, 
more stable and equally as accurate as the other 
methods. The DDM has been combined with the 
LSODE integrator by Dunker 43 into a very efficient 
code. CHEMDDM, which solves both the chemical 
modeling and sensitivity-analysis ODEs for constant 
temperature problems only. Radharkrishnan 44 has 
developed a nonisothermal version of the DDM 
which uses sensitivity analysis subroutines adapted 
from CHEMDDM. These are combined with a 
modified version of the LSODE code in the LSENS 
code mentioned previously. 

The LSENS code 

Description of the code. The LSENS code has been 
designed to model a variety of chemical kinetics 
problems. These include (1) static reacting system, (2) 
steady, one-dimensional reacting flow, (3) shock in- 
itiated chemical reaction, and (4) the perfectly stirred 
reactor. In addition equilibrium chemical compu- 
tations can be performed for four different assigned 
states. Any reaction problem may be adiabatic or 
have a prescribed heat-transfer rate profile. In ad- 
dition, for static or flow reactions, a temperature 
profile may be assigned. In the case of static problems 
either density is held constant or pressure is assigned 
as a function of time. For a flow problem either 
pressure or area may be assigned as a function of time 
or distance. The code is efficient, accurate, and user 
friendly, with several convenience features. 

During a static reaction sensitivity coefficients can 
be computed for all dependent variables and their 
time derivatives with respect to rate coefficient par- 
ameters as well as the initial values of the dependent 
variables. When sensitivity-analysis computations are 
asked for, the code uses the DDM to solve the 
kinetics ODEs, Eq. (15), and the sensitivity coefficient 
ODEs, Eq. (18), sequentially, as mentioned above. It 
should be noted that the coefficients df/dy, in Eq. ( 1 8) 
are exactly the Jacobian elements computed for the 
solution of Eq. (15) by the predictor-corrector tech- 
nique of LSODE. After the kinetics solution is ad- 
vanced from t„ _ [to t„, these results are used with the 
current Jacobian elements to easily solve the sensi- 
tivity coefficient equations from t n _ , to t n . The latter 
solution requires no predictor-corrector procedure. 


Both the kinetics and sensitivity coefficient sol- 
utions of LSENS have been tested by reproducing 
published results from other available codes. 10 Sensi- 
tivity coefficients computed by LSENS were checked 
for several cases by also using the brute-force method 
to calculate the {S,j\. In all comparisons the results 
from LSENS agreed well with the published data and 
with the brute-force results. The computational work 
for LSENS was less than or equal to that needed by 
the other codes for all tests against literature data. 

Application to mechanism development. Research 
has been ongoing for many years to elucidate the 
oxidation mechanism of fuels like hydrogen (H : ) and 
simple straight-chain and cyclic hydrocarbons 
(C,H, ). That this task is very difficult is evidenced by 
the fact that the oxidation mechanism for the simplest 
hydrocarbon, methane (CH 4 ) contains over 125 indi- 
vidual steps. Moreover, the rate coefficient ex- 
pressions for many of the proposed reactions may not 
be known accurately. The approach in this work is to 
propose a set of reactions with either theoretically 
estimated or (often conflicting) literature values for 
the rate parameters, and then attempt to compute 
results obtained in a series of experiments. One needs, 
therefore, not only a chemical kinetics but also a 
sensitivity analysis computational capability to 
efficiently determine the effect of rate-parameter un- 
certainty. Although much information is now known 
about the oxidation of many straight-chain hydrocar- 
bons, 45 the simple aromatic fuels benzene (C 6 H 6 ) and 
toluene (C 7 H 8 ) have not been as well studied until 
recently. These ring structure fuels react differently 
from the straight-chain compounds and now consti- 
tute a relatively high percentage of today’s practical 
fuels. It is, therefore, important to understand the 
oxidation mechanism of these aromatic compounds. 
Recent papers by Bittker 46 and Emdee et aid' have 
presented detailed mechanisms for the oxidation of 
these fuels. We will give a brief summary of the work 
of Bittker on the fuel benzene. First, however, it 
should be pointed out that one does not use a 
completely new mechanism for each fuel studied. 
Because the larger-molecule fuels react to form 
smaller hydrocarbon molecules, the oxidation mech- 
anism of the former is built up from the oxidation 
reactions of the smaller molecules. This is demon- 
strated in Fig. 1 . which shows how toluene reacts to 
form benzene, which then forms many other smaller 
hydrocarbons. The oxidation of all hydrocarbons 
forms molecular hydrogen, so its oxidation is import- 
ant in the reaction of any fuel. 

The benzene mechanism of Bittker employs 29 
reactions of benzene and its fragments and 91 oxi- 
dation reactions of the fuels methane, ethylene 
(C 2 H 4 ), acetylene (C : H : ) and hydrogen to compute 
two different types of experimental data for benzene 
oxidation reported in the literature. The LSENS code 
was used for these kinetics computations and to 
compute sensitivity coefficients for all dependent vari- 
ables with respect to all rate coefficients. This 
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“optimal” mechanism was developed by using the 
results of the sensitivity analysis computations to find 
a small number of reactions whose rate coefficients 
not only had a significant effect on most of the 
experimental results, but also were not accurately 
known. These rate coefficients were then adjusted, 
using sensitivity coefficients as a guide, to give 
the best possible agreement. The signs and magni- 
tudes of the sensitivity coefficients showed when 
further changing of these important rate coefficients 
would not improve the overall agreement. The vari- 
ation of rate parameters stops when any change that 
improves agreement for one variable would degrade 
already good agreement for one or more other vari- 
ables. 

Comparisons of experimental with computed re- 
sults using the optimal mechanism are shown in 
Figs 2-4. In Fig. 2 we present concentration versus 
time profiles for benzene and carbon monoxide (CO) 
which were measured by Lovell et al. 4g in a constant 
pressure, well-mixed flow reactor. Benzene was in- 
jected into a flowing dilute mixture of oxygen in 
nitrogen carrier gas and rapidly formed a homo- 
geneous mixture at a temperature of about HOOK.. 
Velocity profile measurements were used to convert 
any reactor distance from injection point into a 
reaction time, so that the reaction could be modeled 
computationally as a constant-pressure static process. 
Results are shown for two fuel equivalence ratios 
i.e. the ratio of the the actual molar benzene-oxygen 
ratio to the stoichiometric value. The agreement 
between computation and experiment is generally 
good, but is better at the lean 4> value of 0.74 
[Fig. 2(a)] than the rich value of 1.36 [Fig. 2(b)], Also 
the relative positions of the computed and experimen- 



Fig. 1. Interrelation of hydrocarbon oxidation mechanisms. 




Time, msec 

(b) «p = 1.36, T 0 = 1096 K. 

Fig. 2. Benzene and CO concentration versus time in a flow 
reactor. 

tal curves change for the two mixtures. Figures 3 and 
4 show ignition delay times versus reciprocal Kelvin 
temperature for two benzene-oxygen-argon mixtures 



1515 1430 1350 

Temperature, K 

Fig. 3. Benzene-oxygen-argon ignition delay time versus 
reciprocal of temperature, <p = 1.0, dilute mixture 95.6% Ar. 
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Fig. 4. Benzene-oxygen-argon ignition delay time versus 
reciprocal of temperature, (p = 2.0. 


ignited behind a reflected shock wave. The exper- 
iments were reported by Burcat et al* 9 and the delay 
times were based on pressure versus time measure- 
ments. The reaction behind the shock was modeled, 
to a good approximation, as a constant-volume static 
process, and the computed ignition delay times were 
obtained from calculated pressure versus time plots. 
We see that agreement between experiment and com- 
putation is much better for the stoichiometric 
(4> = 1.0) mixture (Fig. 3) than for the rich (tj> = 2.0) 
mixture (Fig. 4). 

It is apparent from these comparisons that 
our mechanism is still incomplete. Additional 
experimental data are needed so that more compari- 
sons between theory and experiment can be done 
to find as yet undiscovered important reactions 
and obtain good agreement with all experimental 
data. This example illustrates the fact that one 
cannot be completely sure that a chemical mechanism 
is correct and complete even if it successfully predicts 
a wide range of experimental data. In fact the 
mechanism may not be unique. The work of Emdee 
et alf 1 presents a benzene submechanism, as part of 
its toluene oxidation mechanism, which differs in 
several respects from that of Bittker. However, their 
mechanism predicts the flow reactor profiles of 
Lovell et alf* equally as well as Bittker’s does. It 
has not been tested as to its ability to reproduce 
the ignition delay time data of Burcat. 49 In summary, 
obtaining knowledge about the mechanism of the 
oxidation of hydrocarbon fuels, or any complex 
chemical process, is a continuous task. As the 


data base of experimental information about the 
reaction increases, the mechanism may have to 
be modified to be consistent with all the available 
data. ' 

SIMPLIFICATION OF CHEMICAL KINETICS 

SOLUTIONS 

There are two basically different approaches to 
the problem of obtaining realistic, but computation- 
ally rapid, combustion kinetics solutions which 
can be used in CFD modeling codes. The first is to 
reduce the size of the reaction mechanism, while 
keeping intact its important features. The speed of a 
chemical kinetics numerical integration increases 
rapidly as the number of reactions and species is 
decreased. The second is a relatively new idea called 
solution mapping, which optimizes the solution 
with respect to a small number of input variables and 
then replaces the differential equation solutions by a 
set of algebraic equations which express the depen- 
dent variables as functions of certain input variables 
of the problem. We conclude this paper with an 
overview of each of these techniques and their current 
status. 

Reduction of chemical mechanisms 

Algebraic simplification. For many years research 
has been directed towards separating the reactions 
and species in a mechanism into groups determined 
by different species reactivities. The solution can be 
simplified by obtaining algebraic relationships involv- 
ing the concentrations of certain subgroups of very 
reactive species and rate coefficients of their reactions. 
This removes some species as unknowns from the 
differential equation set [Eq. (4)] and allows their 
concentrations to be computed algebraically. The 
numerical solution of the reduced set of species ODEs 
is then much easier and much more rapid. The three 
methods that have been used to accomplish this type 
of simplification are (1) the pseudo steady-state ap- 
proximation, (2) the partial equilibrium approxi- 
mation, and (3) the rate-controlled constrained 
equilibrium method. 

The pseudo steady-state approximation (PSSA) is 
based on the assumption that the rates of formation 
of certain very reactive atoms and intermediate 
species called free radicals are very much smaller than 
those of the other species present. This is often true 
because their absolute concentrations are extremely 
small. The net rates of formation of these trace 
species can then be set equal to zero. They are 
assumed to achieve this psuedo steady state in a time 
interval which is short compared to the total reaction 
time. A simple example will illustrate this technique. 
We consider the reacting system (which may be part 
of a larger mechanism) 

k, kx 

A ^B^±D 

k 2 la 
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where the concentration of B is much smaller than 
that of the other two species. The law of mass action 


gives the following species ODEs: 


d c A 


— =k 2 c B -k i c A 

(20) 

= k \ c A + K c D - (k 2 + k 2 )c B 

(21) 

d c d 


c B — k 4 c D . 

(22) 


We now assume dc g /dr = 0 and solve Eq. (21) for c B . 
When the resulting expression is substituted for c B in 
Eqs (20) and (22), the following set of two differential 
equations in the variables c A and c D is obtained: 

< ^-=(k 2 K i -k l )c A + k 2 K 2 c D (23) 

d ^=k } K l c A + (k } K 2 -k i )c D . (24) 

Here A) and K 2 are functions of the rate coefficients 
/c, , k 2 , k 2 and k 4 . In a complex reaction there can be 
several species which satisfy the conditions for a 
pseudo steady state, and the original ODE set can be 
greatly reduced. The major problem with this tech- 
nique is developing a generalized and efficient method 
of finding the steady-state species and performing the 
reduction for any complex system. Significant pro- 
gress has been made in this task, as shown in the work 
of Benson, 50 Bowen et a/., 51 and Dixon-Lewis et air' 2 
The partial-equilibrium approximation (PEA) ap- 
proach also looks for groups of species and reactions 
which can be used in a partial algebraic solution. 
However, this technique searches for fast reactions, 
which are replaced by algebraic relations involving 
the equilibrium constants rather than rate co- 
efficients. In general these fast reactions are easy to 
identify without having detailed physical insight into 
the reaction mechanism. For example, under combus- 
tion conditions the following three hydrogen-oxygen 
reactions are fast enough to be considered in equi- 
librium: 

h + o 2 ^±oh + o 
o + h 2 ^±oh + h 
h 2 + oh-h 2 o + h. 

Equilibrium theory then gives algebraic relationships 
among these species concentrations under the as- 
sumption that the atoms H and O and OH (hydroxyl) 
radical all maintain their equilibrium concentrations. 
As with the PSSA method, systematic methods are 
needed for finding the partial equilibrium groupings 
in any problem. The work of Dixon-Lewis et al . !2 uses 


both the PSSA and the PEA approaches in their 
method of “composite fluxes.” Recent papers by 
Goddard 53 and Rein 54 report two different general 
methods of separating any reacting system into equi- 
librium and nonequilibrium groupings. The Goddard 
formulation implicitly includes the method of com- 
posite fluxes. The method of Rein uses the atom and 
mass conservation conditions that must be satisfied in 
any chemical system to separate the species and 
reactions into linearly independent and dependent 
groupings, and then applies the partial equilibrium 
approach. 

Another approach for simplifying a complex reac- 
tion mechanism is the rate-controlled constrained- 
equilibrium (RCCE) method of Keck and Gillespie. 55 
In this method only the slower rate-controlling reac- 
tions are treated by finite-rate methods. The method 
is described in detail by Keck S6 and is applied to the 
ignition of hydrogen-oxygen mixtures by Law et al. sl 
The first step in this method is to assume that all 
finite-rate reactions are frozen, derive constraints 
using this assumption and compute the equilibrium 
state that satisfies these conditions. Next, the con- 
straints are recalculated using the rate equations of 
the slow reactions. However, the rate-controlling 
reactions have to be determined by trial and error, so 
this method will not necessarily be computationally 
faster than other methods. 

Quasi-global mechanisms. Another way to reduce 
the number of individual molecular reaction steps in 
a complex mechanism is to combine a group of 
molecular reactions into one overall or global reac- 
tion. The global reaction attempts to summarize the 
effect of several reactions which form and destroy 
atoms and reactive radicals, while changing reactant 
species into stable products. It is chosen to accurately 
predict the formulation or destruction rate of a stable 
species without considering the individual molecular 
steps in the process. An example of this reduction 
technique is given in Fig. 5. We show part of a 
detailed mechanism for the oxidation of hydrogen 
and, below it, a single global reaction that represents 

Molecular mechanism (partial) 

H + 0 2 OH + O 

O + H 2 ^ OH + H 

H 2 + OH ^ H 2 0 + H 

H + OH + M ^ H 2 0 + M (M = collision partner) 

Net effect: 2H 2 + 0 2 => 2H 2 0 


Global reaction 
H 2 + 1/2 0 2 => H z O 

Fig. 5. Molecular mechanism and global reaction for oxi- 
dation of hydrogen. 
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the overall process, which is the conversion of 1 mole 
of H 2 and 0.5 mole of O, to 1 mole of water. It is 
interesting to notice that, if we ignore the fact that the 
four molecular reactions are allowed to be reversible, 
the net effect of adding these reactions is the conver- 
sion of 2 moles of H, and 1 mole of 0 2 to 2 moles of 
water. The global reaction (which will always be 
irreversible) will be an adequate representation of the 
kinetics of the combustion if the major interest is the 
rate of fuel consumption or formation of water. 
However, the rate of this reaction is not determined 
by the law of mass action, but by using available 
experimental data to fit the empirical expression 

dr NS 

-5=t(fin^ <25) 

d / ,Vi 

Here, c, is the concentration of any species in the 
mixture (not only reactants), and k(T) is a tempera- 
ture-dependent rate coefficient in the form of Eq. (2). 
The constants {a,} and the three rate parameters are 
determined by computer trials and fitting of exper- 
imental data. The work of Varma et al . 9 uses the 
global reaction in Fig. 5 to study premixed hydro- 
gen-air laminar flames. 

The rate of a global reaction can depend on the 
concentrations of species which do not participate in 
it, because the reaction lumps together the effects of 
several species, even if they all do not participate 
explicitly in the reaction. The obvious drawback of 
the global-reaction approach is that a rate expression 
determined for one set of experimental conditions 
may not be valid for a different set of experimental 
conditions. Also, for a complex system it is often a 
tedious process to determine how to group the mol- 
ecular reactions and replace some or all of them by 
a small number of global reactions. Nevertheless, this 
technique has been used for many years with some 
success. For example, Peters and Kee 58 and Peters 
and Williams 59 have developed simplified three- and 
four-step mechanisms for the oxidation of methane, 
which have been applied to both premixed and 
laminar diffusion flames. The work of Lam and 
Goussis 60 uses the Computational Singular Pertur- 
bation Technique to systematically divide a complex 
mechanism into reaction groups, each with a single 
characteristic time scale. These groups are used to 
develop a reduced and simplified set of global reac- 
tions. Different mechanisms can be obtained, depend- 
ing on which variables are of interest in a given 
problem. In other work global mechanisms have been 
developed for the combustion of methane 7 and 
propane 61 over a range of mixture equivalence ratio. 
Hautman et al propose a general hydrocarbon 
oxidation global mechanism which was verified over 
a wide range of experimental conditions. In the recent 
work of Kundu and Deur 63 a quasi-global mechanism 
which models nitrogen oxide formation in hydro- 
gen-air combustion was developed. Their reduced 
mechanism, containing both global and molecular 



.4 .5 .6 .7 .8 .9 1.0 

Equivalence ratio 

Fig. 6. NO, emissions comparison as a function of equival- 
ence ratio. 


reactions, successfully computes the results of exper- 
imental nitric oxide emissions measurements of An- 
derson 64 in a laboratory combustor over a wide range 
of fuel equivalence ratio. For these limited exper- 
imental conditions the quasi-global mechanism was 
more successful in modeling the production of nitric 
oxide than was a much larger detailed molecular 
chemical mechanism. These comparisons are shown 
in Fig. 6, which is taken from Kundu and Deur. 

Summary. Many partially successful reduced 
mechanisms, employing all of the simplification tech- 
niques that have been discussed here, have been 
developed. However, none as yet combines good 
predictive capability with significantly increased com- 
putational speed for a very wide range of experimen- 
tal conditions, and research continues in this field. 
For a very good recent review of combustion model- 
ing and the use of simplified reaction mechanisms, the 
paper by Dryer 65 is recommended. 

Solution mapping 

An entirely different approach to the rapid compu- 
tation of the kinetics of a complex reaction is the idea 
of solution mapping, first suggested by Frenklach. 8 - 66 
In chemical modeling, the differential equations de- 
scribing the system generate solutions (concen- 
trations, temperature) which are model responses to 
the input parameters of the problem (initial con- 
ditions, rate coefficients), which are called the model 
variables. The idea of solution mapping is to first 
obtain solutions of the ODEs (responses) for a range 
of variation of the input variables. These solutions 
are called computer experiments, and their number 
does not have to be large if the proper combinations 
of model variables are selected. The relationship 
between the model responses and the model variables 
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is then analysed systematically by assuming the re- 
sponses are primarily dependent on a limited number 
of model variables, which are called active par- 
ameters. The latter can be found by using the method 
of local sensitivity analysis described earlier in this 
paper. The model responses are expressed as alge- 
braic functions (often polynomials) of the active 
parameters. These functions are called response sur- 
faces, and they are a reduced model of the kinetic 
process, giving essential dependencies between the 
variables of the problem. It should be pointed out 
that a model response can be a quantity derived from 
the temperature and concentration versus time sol- 
utions of Eqs (4) and (14). An example of this is 
experimental ignition delay time, which can be 
defined by several different criteria. These include a 
fixed pressure rise, 49 or temperature rise, or the time 
needed for a radical species concentration (e.g. hy- 
droxyl, OH), to reach a prescribed value. 

A detailed description of the solution mapping 
technique is given by Frenklach et al. ,:1 who also 
describe its application to the kinetics of the combus- 
tion of methane. In this work, response surfaces were 
computed for several ignition delay times and one 
radical concentration. Although their detailed 
methane oxidation mechanism contains over 140 
reactions, a systematic sensitivity analysis identified 
only 13 active parameters (i.e. reaction rate co- 
efficients) which controlled the computed results. The 
resulting polynomials were used to optimize the 
chemical mechanism by minimizing a weighted least- 
squares function of the differences between computed 
and experimental results. The authors point out the 
significant savings in computer time achieved by 
substituting evaluations of polynomials for repeated 
solutions of sets of differential equations when mini- 
mizing the least-squares function. 

This procedure has also been applied to the sim- 
plification of photochemical air pollution models, 
which require the solution of large systems of stiff 
ODEs (Marsden et al. 6> ). Algebraic models for ozone 
formation were obtained which computed ozone con- 
centration profiles in excellent agreement with the 
results of complete ODE model solutions. The ap- 
proach used in this work could be applied to a 
combustion flow problem by generating polynomial 
functions for temperature and a limited number of 
species concentrations (e.g. fuel, carbon monoxide, 
nitric oxide) and building these functions into a 
subroutine of a CFD code. Repeated calling of this 
subroutine would clearly be less costly than repeated 
calling of an ODE solver subprogram. 

CONCLUDING REMARKS 

In this paper we have discussed the need for 
incorporating realistic heat release and species for- 
mation rate computations into CFD reacting flow 
codes which model practical combustion systems. 
The general methods by which we theoretically pre- 


dict the effects of complex chemical reaction in the 
oxidation and combustion of fuels were described. 
Because the detailed solutions of large chemical sys- 
tems are computationally expensive, they must be 
simplified in order to be used as part of a large 
reactive-flow modeling code. Several methods of sim- 
plification were described. They include (1) using 
pseudo steady-state or partial equilibrium approxi- 
mations to reduce the number of differential 
equations to be solved and substitute algebraic re- 
lations, (2) forming a global or quasi-global reduced 
mechanism which preserves the important features of 
the detailed mechanism, and (3) using the technique 
of solution mapping to transform the results of a 
series of detailed computations into algebraic 
equations for dependent variables of interest. These 
equations are functions of several important input 
“active parameters,” which are found by a careful 
sensitivity analysis of the detailed solutions. Each of 
these simplification approaches has shown good re- 
sults when applied to certain selected systems. How- 
ever, further research is needed on all these 
techniques to develop an efficient and very general 
method for chemical kinetic computations in CFD 
reactive flow codes. 

Acknowledgement — The author thanks Dr K. Radhakrish- 
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